function Lrk4_sin

%  Use RK4 (in t) and centered diff (in x) to solve the heat equation for various M values
%       diff(u,x,x) = diff(u,t)   for xL < x < xr, 0 < t < tmax
%  where
%      u = 0  at x=xL,xR  and  u = sin(2*pi*x) at t = 0

% clear all previous variables and plots
clear *
clf

% set parameters
N=18;
M=5;
tmax=0.1;
xL=0;
xR=1;

% pick time points (by picking the index)
itotal=3;
it(1)=2;
it(2)=(M+1)/2;
it(3)=(M+1);

% generate the points along the x-axis, x(1)=xL and x(N+2)=xR
x=linspace(xL,xR,N+2);
h=x(2)-x(1);

% calculate numerical solution
%set(gcf,'Position', [662 315 560 725]);
plotsize(560,725)
for im=1:3

	% generate the points along the t-axis, t(1)=0 and t(M+1)=tmax
	t=linspace(0,tmax,M+1);
	k=t(2)-t(1);
	lamda=k/h^2;

	fprintf('\n    Lamda = %5.2e\n\n',lamda)
	if im==1
		ue=rk4(x,t,N+2,M+1,h,k);
		tt(1)=t(it(1)); tt(2)=t(it(2)); tt(3)=t(it(3));
	elseif im==2
		uee=rk4(x,t,N+2,M+1,h,k);
	else im==3
		ueee=rk4(x,t,N+2,M+1,h,k);
	end;
	M=2*M;
end;

% plot results
xx=linspace(xL,xR,100);
for itt=1:itotal

	% plot numerical solutions
	subplot(3,1,4-itt)
	hold on
	plot(x,ue(:,it(itt)),'-sr')
	plot(x,uee(:,2*it(itt)-1),'-ob')
	plot(x,ueee(:,4*it(itt)-3),'--','Color',[0.5 0 0.5],'Linewidth',1)
	%plot(x,ueee(:,4*it(itt)-3),'--b','Linewidth',1)
	
	% plot exact solution
	u=exp(-4*pi*pi*tt(itt))*sin(2*pi*xx);
	plot(xx,u,'-k')
	
	% define axes used in plot
	xlabel('x-axis','FontSize',14,'FontWeight','bold')
	ylabel('Solution','FontSize',14,'FontWeight','bold')
	
	% have MATLAB use certain plot options (all are optional)
	set(gca,'FontSize',14); 

	box on
	say=['Time = ', num2str(tt(itt))];
	if itt==1
		yt=0.39;
		axis([0 1 -0.48 0.48]);
		set(gca,'ytick',[-0.48 -0.24 0 0.24 0.48]);
		legend(' M = 5',' M = 10',' M = 20',' Exact',3);
		%legend(' M = 5',' M = 20',' Exact',3);
		set(findobj(gcf,'tag','legend'),'FontSize',12,'FontWeight','bold'); 
	elseif itt==2
		yt=0.17;
		axis([0 1 -0.22 0.22]);
		set(gca,'ytick',[-0.22 -0.11 0 0.11 0.22]);
	else
		yt=0.8*1e20;
		%axis([0 1 -0.02 0.02]);
		%set(gca,'ytick',[-60  -30 0 30 60]);
	end
	text(0.75,yt,say,'FontSize',14,'FontWeight','bold')
	hold off

end;
say=['Heat Equation: exact vs L-RK4 when u(x,0)=sin(2\pix)'];
title(say,'FontSize',14,'FontWeight','bold')


% rk4 method
function UC=rk4(x,t,N,M,h,k)
UC=zeros(N,M);
for i=1:N
	UC(i,1)=g(x(i));
end;
N2=N-2;
alpha=1/(h*h);
A=diag(-2*alpha*ones(N2,1))+diag(alpha*ones(N2-1,1),1)+diag(alpha*ones(N2-1,1),-1);
y=zeros(N2,1);
for i=1:N2
	y(i)=g(x(i+1));
end;
size(y);
size(UC(:,1));
for j=2:M
	k1=k*A*y;
	k2=k*A*(y+0.5*k1);
	k3=k*A*(y+0.5*k2);
	k4=k*A*(y+k3);
	y=y+(k1+2*k2+2*k3+k4)/6;
	UC(:,j) = [0; y; 0];
end;

% subfunction g(x)
function q=g(x)
q=sin(2*pi*x);

% subfunction plotsize
function plotsize(width,height)
siz=get(0,'ScreenSize');
bottom=max(siz(4)-height-95,1);
set(gcf,'Position', [2 bottom width height]);